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1.  INTRODUCTION 

The  statistical  properties  of  estimators  of  signal-to-noise  ratio  of  acoustic 
tones  are  needed  in  ongoing  studies  of  acoustic  signal  processing  (Task  No.DST 
79/069  -  Signal  Processing  for  Underwater  Detection)  in  areas  as  widely  separated 
as  measurement  of  acoustic  transmission  loss  and  signal  detection  theory.  In  a 
previous  paper  by  the  author (ref.  1)  an  expression  was  derived  for  the  probabil¬ 
ity  density  function  of  an  estimate  of  the  signal-to-noise  ratio  of  a  sine  wave 
in  Gaussian  no^t-e.  The  estimate  was  obtained  from  the  power  spectrum  of  a 
sampled  data  sequence.  The  analytic  results  in  reference  1  have  been  used  to 
obtain  functional  forms  for  the  coefficient  of  variation  of  signal-to-noise  ratio 
estimates  derived  from  data  obtained  from  several  sea-going  experiments.  The 
agreement  between  analysis  and  experiment  was  found  to  be  excellent. 

Reference  1  did  not  address  the  situation  where  the  spectrum  used  for  estimation 
was  obtained  as  the  average  of  a  number  of  spectra  derived  from  consecutive 
sequences  of  data  samples.  This  present  paper  gives  a  solution  to  this  problem 
by  using  a  characteristic  function  approach  with  the  aid  of  some  integrals  presented 
in  Appendix  I.  The  problem  is  to  find  the  distribution  of  z 


is  the  power  in  the  j-th  noise  bin  for  the  i-th  spectral  estimate. 

The  noise  bin  must  not  include  the  signal  bin  or  contain  extraneous  signal 
residues. 

N  is  the  number  of  spectra  from  which  an  average  spectrum  is  obtained, 
p  is  the  number  of  noise  bins  in  the  average  spectrum  used  to  estimate  noise 
power 

i  is  the  bin  width  in  Hertz. 

Equation  (1)  can  be  rewritten  as 
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Before  proceeding  with  the  development  of  the  solution,  some  results  of  a  standard 
nature  are  summarised  from  reference  2  in  order  to  aid  in  reading  this  paper,  viz: 

(1)  The  characteristic  function  of  a  random  variable  x  with  probability 
density  function  f(x)  is  given  by 

OO 

^  (t)  J  dx  exp(jtx)  f(x)  (3) 

.OO 

(2)  The  characteristic  function  of  the  sum  S  of  N  random  variables  whose 
individual  characteristic  functions  are  <t>( t)  is 


0s(t)  =  {«( t)}N 


(3)  If  0(t)  is  absolutely  integrable  over  the  range  (-°°,  *9  then 


f(x) 


OO 

dt  exp(-jtx)  0(t) 

.OO 


(4) 


(5) 


2.  THE  CHARACTERISTIC  FUNCTION  OF  SIGNAL  PLUS  NOISE  POWER 

For  a  single  spectral  estimate  the  probability  density  function  of  the  signal  + 
noise  power  is  given  by  equation  (7)  of  reference  1.  Note  that  the  nomenclature 
is  the  same,  ie: 


f(x)  =  A2  exp  l  -A*  (x  +  K2)]  IQ(2A2  K/O,  x  >  0  (6) 


=  0,  x  <  0 

where  Iq(-)  is  the  modified  Bessel  function  of  zero  order.  The  characteristic 
function  of  this  density  function  is  then,  after  a  simple  rearrangement,  given  by 

OO 

*(t)  =  A2  exp(-A2  K2)  J  dx  exp  [-(A2  -  jt)  x  j  IQ(2A2  K  v6T )  (7) 

0 

By  making  the  elementary  substitutions 

o  =  A2  exp(-A2  K2) 

P  =  A2  -  jt 

y  =  2j  A2  K 

JU  =  V$T 


(8) 
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and  using  the  relation  IQ(fO  =  JQ(jM)  the  characteristic  function  can  be  written 
as 


<t>(t)  =  2a  [  n  exp(-/3  M2 )  Jq(7*0 


(9) 


Applying  the  formula  1.1  from  the  Appendix  I  gives 


0(t)  =  f  exp  (-^) 


CIO) 


ie 


0(t) 


,2  1  A 

A  )  F 


-i-  it! 

I-P-} 


exp ( - A 2  X  )  exp 


h* 


K’lfiJi 

A  +  t2 


(11) 


3.  THE  PROBABILITY  DENSITY  FUNCTION  OF  SIGNAL  PLUS  NOISE  POWER 

The  density  function  derived  in  this  section  is  for  the  summation  in  the  numerator 
of  equation  (2). 

For  the  sum  of  N  terms  each  distributed  as  in  equation  (6) 


0sum  =  =  [ff  eXP  [w  ] 

....  from  equation  (10) 

Hence  the  probability  density  function  for  the  sum  is 


f(x) 


1_ 
27 T 


exp(-jt  x) 


exp 


(12) 


(13) 


....  from  equation  (5) 

As  0  is  a  function  of  t,  viz  P  =  A2  -  jt,  the  integral  can  be  rewritten  in  terms 
of  P  to  give,  after  some  rearrangement: 

A2  -  j°° 

f(x)  =  exp (-A2  x)  /  d/5  i  exp(x/J  - -j£)  (14) 

2  it  J  p™  *** 

A2  +  j«* 

Using  the  formula  1.2  in  the  Appendix  I  gives,  after  some  algebra: 

f(x)  =  “^1  V^)(N"1)/2  exP'A2(x+NK^ 

_ _ _ _ _ _ _ _ _ A 
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The  Bessel  function  of  the  first  kind  of  order  n  is  related  to  the  modified 
Bessel  function  of  order  n  by 


t 

E. 

\ 

L 

i 


;  i 

i 

i 


In(x)  =  (-j)n  Jn(jx) 


Hence 


fw  - 


exp  |  -A2  (x  +  NK2  )}  I  j  (2A2  K  \^x)  (15) 


4.  THE  PROBABILITY  DENSITY  FUNCTION  OF  THE  NOISE  POWER 

From  reference  1,  equation  (8),  each  term  y  .  in  the  denominator  of  equation  (2) 
is  distributed  as  X2  2  • 

Hence,  using  the  well-known  addition  theorem(ref . 2)  for  X2  variates,  the  double 
summation  term  in  equation  (2)  is  distributed  as  X2  2p^  1 


if 


M  = 


EE 

j=l  i=l 


13 


then 


a2Pn i2Pni yPN_1 _  exp (-A2  y) 

2pN.  r(pN) 

(16) 

A2PNyPN  I.exp(-A2y) 

P(pN) 


5.  THE  PROBABILITY  DENSITY  FUNCTION  FOR  THE  SIGNAL-TO-NOISE 

RATIO  ESTIMATOR 


the  ratio  of  the  summation  terms  in 
applying  the  classical  relation  given  in 
ratio  of  two  independent  positive  random 

OO 


The  probability  density  function  for 
equation  (2)  can  now  be  obtained  by 
Parzenfref . 2)  for  the  density  of  the 
variables,  viz: 

fx/YM  ‘  j 

0 


fx(yx)fy(x)  dx  (17) 
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le 


t'x/y(y)  =  /  tlx  X  A' 


*2  m 


(N-D/2 


Nfc 


exp(-A2 (yx+NK  )) 


IN1(2A2K  v^xy)  xpN_1  expC-A^x) 

T(pN) 


(18) 


Gathering  terms  independent  of  x  gives 


T(pN)  '•NK' 

The  integral  rearranges  to 


A2pN+2  (•  V  ^  ^N-l)/2  exp(_NA2K2j 


(19) 


OO 

J  dx  xpN+^N  1^2  exp(-A2x(y+l))  1^  ^  C2A2  l&yfjxy) 


(20) 


As  stated  previously  a  simple  relation  holds  between  the  modified  Bessel  functions 
and  the  Bessel  functions  of  the  first  kind,  ie: 


Ik(x)  =  j"k  Jk(jx) 


(21) 


Hence  the  integral  becomes 

OO 

j  dx  xpN+^N  1^2  exp(-A2x(y+l))  j  N+1  JN1(2A2Kj  vffiyx) 
0 

With  the  substitution  P2  =  x  the  integral  becomes 

OO 

2j'N+1  J  dp  p2pN+N  exp(-A2  (y+1)  P2 )  (2A2  Kjv/fty  M) 

0 

=  2j'N+1  j  d  P  P0'1  exp(-02  P2)  Jn  l(7  P) 

0 


(22) 


(23) 


where  the  substitutions 
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a  =  2pN+N+l 
U  =  ACy+l)^ 

7  =  2A2  K  jV$y 


are  distinct  from  the  use  previously  made  in  manipulating  equation  (7). 

A  simple  application  of  the  formula  given  by  equation  1.1  in  the  Appendix  I 
enables  the  integral  (22)  to  be  evaluated  to  be  evaluated  to  give 


„ . -N+l  „,a+n-l.  N-l  ..,a+N-l  „ 

-J  U — — )  7  M(  ^  *4^') 


T(N)  2N  |3a+N'1 


Substitution  for  a,  /?,  7  gives,  after  some  rearrangement: 


(24) 


r (pN+N)  1^  1  N^/2  y  fN-1)/2  M(pN+N,N,  NA2K2  (^j)) 
r(N)  A2pN+2  (y+l)pN+N 

Forming  the  product  with  (19)  gives,  after  some  algebra: 

f  (y)  =  exp(-NA2K2)_y^^  M  (pN+N ,  N ,  NA2  K2  (n^j)) 

X/Y  0(pN,N)  (y+l)pN  N  7 


(25) 


(26) 


where  0(...)  is  the  bivariate  P -function. 

The  distribution  of  z  as  defined  by  equation  (2)  then  follows  from 

r  r,  (  rz+l\ 


fZ(2)  -  tfv(^) 

where  z  =  ^  y  -  —  from  equation  (2). 


(27) 


The  above  is  simply  derived  from  the  classic  relation(ref . 2) 

f  (yl  =  i  f  (Zll2-) 

raX+Biyj  a  X  1  a  ’ 


(28) 


>2„2 


le 


fz(2)  - 


ppN  e'NA  K  (rz-H)N~l  M(pN+N,N,NA 
P (pN,N)  (rz+p+l)pN+N 


(29) 


A  simple  consideration  of  equation  (2)  indicates  that  z  will  only  take  values  in 
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6.  THE  MEAN  AND  VARIANCE  OF  THE  SIGNAL-TO-NOISE  RATIO  ESTIMATOR 

The  mean  and  variance  of  z  follow  readily  from  the  moments  of  y,  where  y  is  the 
ratio  of  the  summation  terms  appearing  in  equation  (2).  The  moments  of  y  can  be 
derived  from  the  density  equation  (26),  ie: 


E(  ynl  = 


exp  (-NA2  K2  ) 
0(pN,N) 


n+CN-1) 


dy  (y+1)pN+N  MCpN+N.N.NA2^^)) 


A  series  expansion  for  M written  in  the  form 


u,  ,  ,  T(b)  \  T(a+r)  zr 
M(a,b, z)  -  Tfl)  /  f (b+r)  r! 


can  be  substituted  in  the  above  equation  to  give,  after  reversing  the  order  of 
integration  and  summation  and  using  the  T-function  form  of  the  0-function: 


E(pN+N+r)  (NA2  K2  )  r  f 
r(N+r)  r!  I 


T(pN) 


n+N-l+r 

_y _ 

dy  (y+i)pN+N+r 


As  shown  in  the  derivation  1.3  in  the  appendix  the  integral  is  simply  the 
0-function 

0 ( (n+N+r) , (pN-n) ) 

This  0-function  can  be  written  in  terms  of  appropriate  T-functions  and  after  a 
rearrangement  of  terms  the  n-th  moment  of  y  can  be  written  as 

Ely11]  =  exp(-NA2K2)  M(N+n,N,NA2  K2  )  (33 

Hence 

Ely]  =  -  exP(-NA2K2)  M(N+!  ,N,NA2  K2 )  (34 

'  pN-1 

A  simple  application  of  the  formulae  13.4.1  and  13.6.12  in  Abromawitz  and  Stegun 
(ref. 3)  gives 

=  n(UA2k2) 


A  similar  application  of  the  same  formulae  enables  the  second  moment  of  y  to  be 
written 


Ely2! 


(pN-1) (pN-2) 


{  (NA2  K2  +N+2)  (1+A2  K2  )  -lj 


] 
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From  equation  (2) 


and 


El  2]  =  £  Ely]  -  1 

r  r 


var(z) 


(37) 


(38) 


Hence  precise  expressions  can  be  written  for  the  mean  and  variance  of  z,  ie: 


E[  2] 


pN(l+A2  K2  )  1 

r  (pN-1)  ~  r 


(39) 


and 


var(z) 


+  2NfpN+N-l)A2  K2  +  (pN2  +N2  -N) 
(pN-1)2  (pN-2) 


(40) 


And,  the  coefficient  of  variation  defined  as 


where  a  is  the  standard  deviation  can  be  written  down  immediately. 
Substituting  N=1  in  equations  (39)  and  (40)  gives 


(41) 


E[z]  =  -  I 

r(p-l)  r 

varU)  -  41  AV  *  *  p)% 

7l  P-2  J 

and  after  some  algebra 

p  ;  AV+  2pA2  K2  ♦ 

C  "  pA2K2  +  1  L  p-2  J 


These  are  precisely  the  respective  equations  developed  in  reference  1  and  hence 
equations  (42),  (43)  and  (44)  supply  a  partial  validation  of  the  analysis 
contained  in  this  memorandum. 


(42) 

(43) 
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7.  AN  ALTERNATIVE  DERIVATION 


On  realising  (note  the  acknowledgement  at  the  end  of  this  paper)  that  the  quotient 
in  equation  (1)  is  the  ratio  of  a  non-central  chi-squared  variate  (with  2N  degrees 
of  freedom)  to  a  chi-squared  variate  with  2pN  degrees  of  freedom  ie  the  quotient 
has  the  distribution  of  a  non-central  F  distribution,  it  is  possible  to  derive 
the  formulas  equations  (29),  (38)  and  (40)  as  particular  cases  of  a  more  general 
formula.  The  non-central  F  with  degrees  of  freedom  in  the  numerator  and 
v2  degrees  of  freedom  in  the  denominator  has  a  probability  density  function 


fY(y) 


exp(-\/2)  l  Vx  ‘'i  /2  Ui /2  -  1 
r(«'a/2)  [i'jJ 


X  1  zT  <>*  +  ^/2  x 

L  "j  J 


T  ILLl 

)  TtyTTTTj)  jt 


i) 


x 


(45) 


where  X  is  the  non-centrality  parameter.  This  equation  can  obviously  be  written 
as 


fY(y) 


ffP.C-X/21 

f^a/2) 


vx/2 


V\/2 


1 


xj  1  +JLL  2 1  ~(yi  +  Vl)/2  rCCy»  v2 )/2) 

l  j  r^j/2) 


fvt  +  V, 


zXp,  \ 
2^2  +  2v\  z ) 


On  making  the  substitutions 


(46) 


vx  =  2N,  v2  =  2pN,  X  =  2NA2  K2 


and  carrying  out  some  elementary  algebra 


fy(y)  =  y 


N  -  1 


ppN  M(pN 


/J(pN,N)  (y  ♦  p) 


♦  N,N, 

pr- 


~nr 


yNA2  K2 

p*y 


(47) 


In  this  alternative  derivation  z  =  1_  y  -  1_ 

r  r 
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Hence 


fz(z)  =  rfy(rz  +  1) 


as  a  result  of  applying  equation  (2) 
dN 


le 


fz(z)  - 


rP 


exp(-NA2  K2 )  (rz  +  1)N  ‘  1  M(pN+N,N,NA2 

j3  (pN,N)  (rz  +  p  +  1)PN  +  N  ~ 


which  agrees  with  equation  (29). 

The  nth  order  moment  of  the  non-central  F- distribution  to  be  derived  from 
equation  (45)  involves  evaluating  an  integral 


/ 


,  n  v2/2~l 
dy  y  v 


vi  ^(^i  +  ^)/2 


(1  +  ~a~  y) 


(l  +  £  y)j 


(48) 


The  integral  can  be  written  as 


OO 

r 


Ci 


dz  - - 

j  (1  +  az) 

0 


where 


(49) 


Ci  =  n  +  j  +  j  -  1,  C2 


+  j,  a  =  -f- 


On  making  the  substitution  az  =  ft  the  integral  becomes 


Ci  +  1 


.Ci 


(1  +  M)Cj 


0(c*  -  c>  -1*  c>  +  *) 

"l  *  1 


(50) 


On  substituting  this  integral  in  the  appropriate  expression  for  E[ y11]  and  carrying 
out  some  algebra 

Ely-1  .  «P(-V2>  &)  *  >4  .  [r(if  •  n) 

v  c  i>2)  r( »,)  J 


i 


11 
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Using  the  substitutions  =  2N,  *  2pN,  X  *  2NA2  KJ  for  values  of  n  of  1 

and  2  in  this  expression  then  enables  El  z]  and  E(  z2 1  to  be  derived  from 


z 


1_ 

r 


y 


i_ 

r 


Although  not  carried  out  in  detail  the  results  are  easily  shown  to  be  identical 
to  equations  (39)  and  (40). 


8.  SOME  PRACTICAL  RESULTS 

In  order  to  give  some  idea  of  how  to  use  the  formulae,  two  steps  were  taken. 

First,  a  program  (program  A  in  Appendix  II)  was  developed  to  examine  the  variation 
in  the  form  of  the  density  function  as  the  number  of  spectra  used  to  obtain  an 
average  spectrum  is  increased.  Results  are  shown  in  figure  1  for  a  moderately 
high  input  signal-to-noise  ratio  of  S  dB,  using  10  noise  bins  and  a  bin  width  of 
1  Hz,  ie  p=10,  r=l .  It  is  surprising  that  even  at  this  signal  level  the  spread 
in  a  number  of  estimates  can  be  high.  This  spread  should  be  compared  with  the 
results  in  figure  2  for  a  much  lower  input  signal-to-noise  ratio  of  -S  dB.  The 
comparison  indicates  a  necessity  for  a  further  detailed  examination. 

The  second  step  taken  was  to  develop  a  program  (program  B  in  the  appendix)  to 
examine  the  variation  with  signal-to-noise  ratio  of  the  mean,  standard  deviation, 
and  coefficient  of  variation  of  the  estimate  for  a  specified  number  of  bins  used 
to  estimate  noise  power.  Typical  results  for  the  mean  and  coefficient  of  vari¬ 
ation  are  shown  in  figures  3  and  4  respectively,  using  p*10  and  r«l  as  in  figure  1. 
In  figure  3  the  bias  (due  to  the  assymetry  of  the  distribution)  is  clear  at  even 
high  signal  levels.  At  low  signal  levels  it  is  obvious  that  it  is  necessary  to 
average  a  large  number  of  spectra  to  minimise  the  bias  in  the  mean.  In  figure  4 
the  behaviour  in  the  coefficient  of  variation  indicates  the  large  number  of  spectra 
to  be  analysed  to  keep  the  frequency  of  occurrence  of  negative  estimates  of 
signal-to-noise  ratio  to  a  minimum.  This  need  is  of  course  going  to  conflict 
with  the  non-stationarity  of  the  statistics  of  any  physical  medium  in  which  an 
experiment  is  taking  place  such  as  when  measuring  acoustic  signal  transmission 
properties  in  the  ocean. 


9.  CONCLUSIONS 

This  paper  presents  a  rigorous  mathematical  analysis  and  a  brief  look  at  some 
experimental  implications  of  a  definition  and  associated  measurement  technique 
of  the  signal-to-noise  ratio  of  a  sine  wave  in  white  noise.  The  computer 
programs  developed  as  a  result  of  this  analysis  can  be  used  to  carry  out  a 
detailed  examination  of  any  proposed  experimental  scenario.  During  the  develop¬ 
ment  of  the  analysis  it  has  appeared  that  the  probability  distribution  function 
for  both  signal  plus  noise  and  signal-to-noise  ratio  are  both  eaenable  to  an 
analysis  that  extends  techniques  reported  in  this  paper  and  it  is  proposed  to 
publish  results  on  this  shortly.  It  is  anticipated  that  this  further  work  will 
throw  some  interesting  light  on  some  of  the  problems  of  detection  theory. 
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APPENDIX  I 

MATHEMATICAL  DETAILS 

1. 1  A  general  integral  formula  first  due  to  Hankel(ref.4) 


dt  t^'1  exp(-p2 12 )  JN(at)  = 


7 - a  M  <£-*-£.  S+l  ,TT) 

r(i»  ♦  i)  2  1  p  M  2  4p 


t  arg  p  I  <  J 


Re(M  ♦  v)  >  0 


If  the  substitutions  a  =  y ,  v  =  o,  m  =  2,  p  =  /jT  are  made  then 

OO 

r  t 

;  dt  t  exp  (-fit2 )  JQ(n)  =  M0,1,  ) 

"o 


1 

2/3 


exP  rf-) 


from  applying  formula  13. 6. 12 (ref. 3) . 

1.2  An  integral  formula(ref .4) ,  generally  attributed  to  Sonine,  for  the  ordinary 
Bessel  function  is 


V 

Chz) 

2irj 


C+j°° 


J 

C+j°° 


exP(t  -  -ft) 


where  the  path  of  integration  is  the  straight  line  R4(t)  *  C  >  0 
If  the  substitutions  »>+l  =  N,  t  =  x/J,  z2  =  Ny2x,  C  =  A2  are  made,  then 

7N-lN(N-l)/2  A2+j°°  ,  ^ 

Jf-i  (7vn^:)  =  ~N  .  (N-l)/2  [  dP^-expC^  -ijf-) 

4  WlX  ‘  " 


Aa-j« 


1.3  0- functions  can  be  written  in  the  integral  form 

OO 

f  yq'1 

*  — — - dy 

;  (y+i)p  q 


0(p.q) 
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Make  the  parameter  substitutions 

q  =  n  +  N  +  r 
p  =  pN  -  n 

Then 

p  +  q  =  pN  +  N  +  r 


Hence 


y.n+N-1+r 

dy 


=  fi  (n+N+r,pN-n) 


as  P  (p,q)  =  0(q»p) 


-  IS 
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APPEND  IX  II 
COMPUTER  PROGRAM 
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Figure  1.  Variation  in  the  probability  density  function  for  the  estimate  of 
signal-to-noise  ratio  from  an  average  of  N  spectra  using  10  noise 
bins  at  an  input  signal-to-noise  ratio  of  5  dB 
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